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Description 

FIELD OF THE INVENTION 

5 [0001] The present invention relates to a method for correcting inhomogeneities in a digital radiation image which 
are inherent to the use of a radiation source, a radiation detector or a radiation image read-out apparatus or the like. 

BACKGROUND OF THE INVENTION 

w [0002] Digital radiography offers the possibility for computer aided diagnosis and quantitative analysis using image 
processing techniques such as segmentation, classification and contrast enhancement. However, computer-based 
image interpretation may be hindered by the presence of the non-uniformities in radiation exposure that are inherent 
to the image formation. 

[0003] In the field of X-ray exposure these non-uniformities can be largely attributed to the Heel effect, nevertheless 

15 other sources of inhomogeneities exist such as recording member non-uniformities or read-out inhomogeneities. 

[0004] Although the intensity inhomogeneities induced by all these factors are smoothly varying as a function of 
location and are easily corrected by the human visual perception system, they complicate the use of automatic process- 
ing techniques because the relative brightness of an object within the image becomes position dependent. The overall 
intensity range is unnecessarily enlarged by the presence of these slowly varying shading components and hence the 

20 dynamic range available to represent diagnostic signal details is reduced. 

[0005] A typical hand radiograph is shown in figure 1 . The background at the left side of the image is clearly brighter 
than at the right side. This phenomenon can be attributed to the so-called Heel effect. Because beam collimator blades 
substantially attenuate the X-ray beam so as to shield irrelevant body parts, the phenomenon is only visible in the 
direct exposure and diagnostic areas and not in the collimation area. The Heel effect can be understood from the 

25 construction of the X-ray tube as schematically depicted in figure 2. Electrons originating from the cathode are attracted 
by the positively charged anode. For better heat dissipation, the anode rotates and is inclined by a small anode angle 
6, which enlarges the area A actua | that is bombarded by electrons while keeping the size of the focal A effl from which 
rays are projected downward to the object, fairly small. As shown in the diagram of figure 3, this design makes the 
length of the path travelled by the X-rays through the anode larger on the anode side of the field T a than on the cathode 

30 side T c . Hence the incident X-ray intensity is smaller at the anode side than at the cathode side of the recording device, 
which explains the inhomogeneity of the background in figure 1 . 

[0006] The Heel effect is one possible cause of intensity inhomogeneities that can be introduced in radiographs, As 
has already been mentioned higher, other causes of non-uniformities might be envisioned such as non-uniform sen- 
sitivity of the recording member, e.g. a photographic film, a photostimulable phosphor screen, a needle phosphor 

35 screen, a direct radiography detector or the like. Still another cause might be non-uniformities of the read-out system 
which is used for reading an image that has been stored in a recording member of the kind described above. 
[0007] Because the image acquisition parameters that affect intensity inhomogeneity vary from image to image (e. 
g. variable positioning of the recording device relative to the X-ray source) and can not be recovered from the acquired 
image at read-out, correction methods based on calibration images or flat field exposure such as the one described in 

40 EP-A-823 691 are not feasible. 

[0008] The method disclosed in EP-A-823 691 comprises the steps of (1) exposing an object to radiation emitted by 
a source of radiation, (2) recording a radiation image of said object on a radiation-sensitive recording member, (3) 
reading the image that has been stored in said recording member and converting the read image into a digital image 
representation, (4) generating a set of correction data and (5) correcting said digital image representation by means 

45 of said set of correction data. The set of correction data is deduced from data corresponding with a uniform, flat field 
exposure of the recording member. The set of correction values represent the deviation of the value that is effectively 
obtained in a pixel of the recording member and a value that would be expected in the case of flat field exposure. 
These correction values associated with each recording member are normally determined once and kept fixed during 
further acquisition cycles. 

so [0009] This type of methods is not applicable for solving problems such as the introduction of inhomogeneities due 
to the Heel effect. 

OBJECTS OF THE INVENTION 

55 [0010] It is an object of the present invention to provide a method to correct a digital image for artefacts such as 
artefacts which are inherent to the use of an X-ray tube, artefacts originating from defects in a radiation detector orthe 
like. 

[001 1] Further objects will become apparent from the description hereafter. 
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SUMMARY OF THE INVENTION 

[0012] The above mentioned objects are realised by a method having the specific features set out in claim 1 . 
[0013] Unlike the state of the art methods, in the present invention a set of correction data is deduced from the actual 
5 image data obtained by exposing an object to radiation and not from an additional image such as an image representing 
a flat field exposure. 

[0014] A radiation image most generally comprises a collimation area which is an area that has been shielded from 
exposure to radiation by shielding elements, a direct exposure area (also called background area) being the area on 
the recording member that has been exposed to direct, unmodulated irradiation and a diagnostic area which is the 
10 area where the radiation image is found of the body that was exposed to radiation. 

[001 5] Because of substantial attenuation of the X-ray beam when passing through the collimator blades, the dynamic 
range of the irradiation has been reduced to an extent so as to make these collimation areas unsuitable for estimation 
of the correction values. 

[0016] Because the causes that introduce inhomogeneities are not present in the collimation area, the collimation 
15 area can be neglected in the method of the present invention. A segmentation algorithm can be used to find the bound- 
aries of the collimation area to exclude these areas from further consideration. An example of such an algorithm has 
been described in European patent applications EP-A-610 605 and EP-A-742 536 (forthecase of a partitioned image), 
these documents are incorporated by reference. 

20 First embodiment 

[001 7] In one embodiment (1 ) a mathematical model representing the phenomenon that induces the inhomogeneities 
is generated. Next (2) the digital image representation is subjected to image segmentation in order to extract data 
representing an estimation of the direct exposure area. Then , (3) parameters of this model are deduced from image 

25 data representing the direct exposure area in the image. Next (4) a bias field is generated on the basis of the deduced 
parameters. Next, (5) a correction by means of said bias field is applied to the image data. Corrected image data are 
then subjected to a stopping criterion. Unless this stopping criterion is met, steps (2) to (6) are repeated. 
[0018] Because the inhomogeneities are only directly measurable in the direct exposure areas of the image, this 
area is preferably first extracted and the parameters of the model are estimated from the data regarding this region 

30 only. A seed fill algorithm can be used to determine the background area. The seed fill algorithm can be started from 
the boundary of the collimation area. 

[0019] Inhomogeneity correction by applying a bias field is performed on the entire image, In the context of the 
present invention the term 'a bias field' is used to denote a low frequency pattern that is superimposed on the average 
image data values in a multiplicative or additive manner. 
35 [0020] Next a new background region is extracted from the corrected image data and the model parameters are re- 
estimated. This sequence is repeated. The method iterates between background segmentation and correction until 
convergence, i.e. until no significant changes in background or parameter estimation occur. 

Second embodiment 

[0021] In a second embodiment according to the present invention a statistical model of the image is first generated 
on the basis of intensity and spatial statistics of image regions in the image. 

[0022] The digital image representation is then subjected to image segmentation in order to extract data constituting 
an estimation of these image regions. The image regions referred to are e.g. direct exposure area, bone image, soft 
45 tissue image etc. 

[0023] Parameters of the statistical model are estimated by means of data pertaining to these image regions. Next, 
a bias field comprising correction data is generated and the entire image is corrected by means of the bias field. The 
result of the previous step is evaluated relative to a stopping criterion. The method steps of segmenting, bias field 
correction and evaluation are repeated until the stopping criterion is met. The stopping criterion is e.g. reached when 
so no significant changes occur in the estimation of the image regions and/or no significant changes occur in the param- 
eters defining the statistical model. 

[0024] In one embodiment the image regions are confined to direct exposure areas. In another embodiment (see 
fourth embodiment) the method steps are applied separately to each of a number of image region classes jointly 
constituting the intensity histogram of the radiation image. 
55 [0025] In one embodiment the statistical model is a Gaussian distribution and the parameters of the statistical model 
are the statistical parameters defining the statistical model such as average value u. of the Gaussian distribution and 
the standard deviation a. 

[0026] The stopping criterion is e.g. reached when no significant changes occur in the estimation of image regions 
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and/or no significant changes occur in the parameters defining the statistical model. 
Third embodiment 

[0027] A third embodiment is based on the observation that the entropy of an image increases if inhomogeneities 
are induced in the image. 

[0028] In a third embodiment of the method of the present invention an information theoretic model of the image 
comprising at least direct exposure areas and diagnostic areas is generated. The model is based on Shannon-Wiener 
entropy increasing when additional intensity value entries are added to the image intensity distribution. The digital 
image representation is subjected to image segmentation in order to extract data representing an estimation of the 
direct exposure areas and the entropy in said model is extracted based on of data pertaining to these areas. Next, a 
bias field is generated and the entire image is corrected by means of the bias field. The result of the previous step is 
evaluated relative to a stopping criterion and the method is repeated until said stopping criterion is met. A stopping 
criterion is e.g. met when the entropy is minimal and no significant changes occur in it. 

Fourth embodiment 

[0029] The fourth embodiment is a more general case of the second embodiment. In this fourth embodiment the 
method steps of the second embodiment are applied separately to each of a number of image region classes jointly 
constituting the intensity histogram of the radiation image. 

[0030] In all of the aforementioned embodiments, the number of iterations may be restricted to one when less pre- 
cision is needed and hence the stopping criterion need not be evaluated. 

[0031] Further advantages and embodiments of the present invention will become apparent from the following de- 
scription [and drawings], 

BRIEF DESCRIPTION OF THE DRAWINGS 

[0032] 

Fig. 1 shows a typical radiographic image of a hand in which the Heel effect is clearly visible on the direct exposure 
area, 

Fig. 2 and 3 are schematic side views of an X-ray tube, 

Fig. 4 is a coordinate system wherein an X-ray originates at position (0,0) and travels along R to the recording 
medium at position (p, DJ. 

DETAILED DESCRIPTION OF THE INVENTION 

First embodiment 

[0033] A mathematical model for the Heel effect can be derived from the simplified one-dimensional model of the 
anode and beam geometry depicted in Fig. 4. In the coordinate system (p,z), with p along the anode-cathode axis and 
z along the vertical direction, the X-rays can be taught off to originate within the anode at point o>(0,0), at a distance 
D ave from the anode surface S. Consider the ray R at an angle e)> from the vertical within the plane (co,S) that hits the 
recording device at point {p,D is ) with D /s the distance between the X-ray source and the recording device and 
tan<|)=j£-. The distance /-traveled by R through the anode is given by 




(1) 



with c,(p R ,z R ) the intersection of ffwith S which can be found by solving the system of equations: 



S: p R -D^ -tan0.Zfl 
R:p R = tan <f>.z R 



(2) 
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[0034] Hence, 




(3) 



[0035] The radiation received on the recording device is 



M{p)=l 0 .e 



(4) 



with n the attenuation coefficient of the anode material and l a the radiation originating at to. 
[0036] Model (4) predicts that the Heel effect behaves exponentially along the anode-cathode axis and assumes 
that it is constant perpendicular to this axis. This is justified by flat field exposure experiments which show that the 
difference in intensity perpendicularto the anode-cathode axis is relatively small compared to the intensity differences 
along the anode-cathode axis. 

Image segmentation: 

[0037] A typical hand radiograph, as shown in Fig. 1 , consists of three regions: collimation area (numeral 1), direct 
exposure area (numeral 2) and diagnostic regions (numeral 3). Because the Heel effect is largely reduced in the col- 
limation area and directly measurable in the direct exposure area only, the image needs to be segmented to fit model 
(4) to the image intensity data. This is obtained by first extracting the collimation area and then searching the direct 
exposure area, the remaining areas being diagnostic regions. 

[0038] The boundaries of the collimation area have been found using the Hough transform, assuming that these are 
rectilinear edges as is the case for the majority of X-ray source-mounted collimation devices. To make this approach 
more robust, the contributions of each image point to the Hough space accumulator are weighted by said point's gradient 
magnitude and, for each point, only lines the direction of which is within 10 degrees from the normal to the local gradient 
direction are considered. The 4 most salient points in Hough space that represent a quadragon with inner angles 
between 80 and 100 degrees are selected as candidate boundaries of the collimation area. Because not all 4 collimation 
shutter blades leave an imprint in the image and hence make the associated boundaries disappear in the image, 
candidate boundaries along which the image intensity differs from the intensity expected for the collimation region are 
rejected. 

[0039] To extract the background region 6, a seed fill algorithm has been used that starts from the boundary of the 
collimation region as determined in the previous step. Appropriate seed points for B are found by considering a small 
band along each of the collimator edges and retaining all pixels whose intensity is smaller than the mean of the band. 
This approach avoids choosing pixels that belong to the diagnostic region as candidate seed pixels. B is then grown 
by considering all neighboring pixels n,-,/=1 ,...,8 of each pixel peB and adding q,to B if the intensity difference between 
p and q, is smaller than some specified threshold. 

Heel effect estimation: 

[0040] To fit the model (4) to the image data N (x,y) the direction y has to be found of the anode-cathode axis and 
the parameters a = [l 0 ,n.,B,D is ,D ave ,p m ] such that the model best fits the image data within the direct exposure area 
extracted thus far. p ro is a parameter introduced to map point to where the X-ray originates to the correct image coor- 
dinates (see Fig. 4). Forthe case whereby the Heel effect is modulated as a one-dimensional phenomenon, the distance 
p m and the angle y are the required parameters to map the coordinate system attached to the X-ray origin w to the 
image plane coordinate system. However, because the anode has the three-dimensional shape of a cone, the Heel 
effect is a three dimensional phenomenon , in that the intensity also slightly is reduced in a direction perpendicular to 
the (p,z) plane. To model the two-dimensional Heel effect in the image plane, a third geometry parameter p ro± is needed. 
Parameters (p^p^y) jointly define a coordinate system translation and rotation from the X-ray origin oo to an image 
plane origin, which is the center of the image e.g. in practice the Heel effect inhomogeneity in said perpendicular 
direction is only small with respect to the Heel effect along the anode-cathode axis. 

[0041] Assuming that y is known, the average image profile P^p) along this direction in the direct exposure region 
6 is given by 
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with x and y the image coordinates as defined in Fig. 3 and <.} the averaging operator. We can then find the 
optimal model parameters a* by fitting the expected profile M (p,a) to the measured profile: 



a" (y) = arg mm||P y [p)-M{p, a| 



[0042] The fitted one-dimensional model M(p,a*{y)) is then back projected perpendicular to the projection axis yto 
obtain a reconstruction R(x,y,y,a\y)) for the whole image: 

20 

R(x,y,y,a *(y)) = /W(x.cosy+y.siny,a *(y)) 

[0043] The direction of the anode-cathode axis y is then determined such that this reconstruction best fits the actual 
25 image data within the direct exposure region using 



r* - arg min || N{x, y )- R{x >y , r ,a{y))\^ 



. N(x,y) 

r = arg mm ^ 1 

r \\R(x,y,r,oc (7)) I 



depending on whether we wish to use additive or multiplicative correction. The estimated Heel effect is R(x,y,f, 
a*(Y)) and the corrected image is respectively 

N(x,y)= N(x,y)-R(x,y,y *« *(y*)) (8) 



"™'-R(x,y,y*,a'(y*)y 

[0044] The optimal parameters a* and y* are found by multidimensional downhill simplex search. It has been noticed 
that the anode-cathode axis in our setup is almost always parallel to the image or collimation edges. This reduces the 
number of orientations which have to be evaluated in (6-7) and reduces computation time. 

[0045] After inhomogeneity correction of the image using (8-9), the direct exposure area Sis updated by thresholding, 
using a threshold derived from the histogram of the corrected image intensities N. Keeping the previously determined 
anode-cathode orientation y, new values for the optimal model parameters a* are determined using (5) taking the newly 
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selected direct exposure region into account. A number of iterations, typically three or four, have been performed 
between background segmentation and Heel effect correction until convergence. 

Second and third embodiment 



[0046] In ideal circumstances, the image formation process of diagnostic digital X-ray images is usually well described 
by a multiplicative model yielding an intensity-uniform image U(x,y): 



where 0(x,y) represents the object in the image. In diagnostic X-ray images, the most important contributing 
process of the object is the linear attenuation of the X-rays by the bone and soft tissue 



u. is the linear attenuation coefficient along the path between the origination X-ray at position to and the recording 
device However, nonuniform illumination / = l(x,y), uneven sensitivity of the recording device and inhomogeneous 
sensitivity of the phosphors for readout, introduce unwanted intensity modulations in the acquired image N(x,y) de- 
scribed by function f 

N{x,y)=f x „ u (U(x,y)) (10) 

[0047] In the second and third embodiment the Heel effect is again examined as a very important source of nonu- 
niform illumination. Reference is made to figures 2 - 4 which aid in explaining this effect. 

[0048] Electrons originating from the cathode are attracted by the positively charged anode. For better heat dissipa- 
tion, the anode rotates and is inclined by a small anode angle 9, which enlarges the area /4 ac([ , a/ that is bombarded by 
electrons while keeping the size of the focal spot A eff , from which rays are projected downward to the object, fairly 
small. As shown in the Fig.3, the design makes the length of the path travelled by the X-rays through the anode larger 
on the anode side of the field (T a ) than on the cathode side (T c ). Hence the incident X-ray intensity is smaller at the 
anode side of the recording device. A simple theoretical model is given by 



with l 0 the radiation originating at w,u, the linear attenuation coefficient of the anode, D ave the average distance 
traveled through the anode by the electrons, D js the distance between the X-ray source and the recording device and 
pthe distance from the recording device to X-ray source projected onto the anode-cathode axis. 
[0049] Although the second and third embodiments are explained with reference to the Heel effect, other source of 
inhomogeneities may be envisaged such as the moulding process of imaging plates and/or the characteristics of the 
read-out system. In some fabrication processes, the concentration of phosphors at the edge of the plate is lower than 
the concentration in the middle of the plate which may result in a non-uniform image. I n read-out apparatuses comprising 
mirror deflection, the displacements of the mirror has to be very accurately implemented to achieve uniform activation 
of the phosphors for read-out. Due to all these factors, it is almost impossible to model the induced inhomogeneities 
correctly and more general image formation models are needed. 

Problem formulation : 

[0050] The image formation process is generally modeled with a function /applied to an ideal intensity-uniform image 
U(x,y), resulting in the acquired image N(x,y). In digital X-ray images, the image degradation process dependency on 
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the intensity values U{x,y) is relatively small compared to position dependent factors. Hence, we can rewrite equation 
(10) as follows 

N(x,y)=f xy (U(x,y)) 

[0051] This equation can be simplified as 

N(x,y)=U(x,y)S^x,y)+S A (x,y) 

where S^x,y) and S A (x,y) represent the multiplicative and additive components of the image degradation proc- 
ess. To remove the image inhomogeneities, a corrected image U\s searched which optimally estimates the true image 
U. If the estimates S A and S M of the actual formation components S A and S M are available, the corrected image 0 is 
given by the inverse of the image formation model 



S„<x,y) 



= N(x,y)S M (x,y)-S A (x,y) 



with S M (x,y)=-^- and SAx,y) = . 

S M (x,y) S M {x,y) 

The problem of correcting the inhomogeneities is thus reformulated as the problem of estimating the additive and 
multiplicative components s,ands„. 

Correction strategy: 

[0052] Finding the optimal parameters of the components s,,and s M involves defining a criterion which has to be 
optimized. In this section, two criterions are defined. 

[0053] One correction strategy (second embodiment of the method according to the present invention) is based on 
the assumption that the intensity values of the direct exposure area (also referred to as background) from the acquired 
image is Gaussian distributed. In ideal circumstances, this assumption is true for the acquired image N(x,y). The like- 
lihood that a pixel u,- of the corrected image belongs to the background is 



4lna 2 



{ 2 a 2 



where u. and a 2 are the true mean and variance of the Gaussian distribution of the background pixels. Given an 
estimate B of the direct exposure area, we seek to maximize the likelihood 
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s which is equivalent to minimizing the log-likelihood 



U* = arg min- Z log e ^(i^a). ( 13 ) 

BJJ ieB 

10 

[0054] Another embodiment (third embodiment of the method of the present invention) is based on the assumption 
that the information content of the acquired image is higher than the information content of the uniform image, due to 
the added complexity of the imposed inhomogeneities: 

l c (N(x,y)) =l c (f Ky U(x,y))) > l c (U(x,y)) 

[0055] The information content l c can be directly expressed by the Shannon-Wiener entropy 



I c (V(x,y» = H(N(x,y))=-Sp( n )lo g . P (») (14) 



where p(n) is the probability that a point in image N{x,y) has intensity value n. The optimal corrected image 0* 
is thus given by 



U* = argmjn H(U(x,y)) ( 15) 

U 



Method 

35 

\0056] Because the Heel effect is almost totally reduced in the collimation area and an estimate of the background 
B is needed to optimize equation (1 3), a segmentation algorithm is presented. 

[0057] In the next, implementation details of the correction models of the second and third embodiment of the method 
according to the present invention are given. 

Image Segmentation 

[0058] The boundaries of the collimation area have been found using the Hough transform, assuming that these are 
rectilinear edges as is the case for all hand radiographs in our database. To make this approach more robust, the 

45 contributions of each image point to the Hough accumulator are weighted by its gradient magnitude and, for each point, 
only the lines whose direction is within 10 degrees of the normal to the local gradient direction are considered. The 4 
most salient points in Hough space that represent a quadragon with inner angles between 80 and 100 degrees are 
selected as candidate boundary of the collimation area. Because not all 4 collimation boundaries are always present 
in the image, candidate boundaries along which the image intensity differ from the expected intensity values for the 

so collimation region, are rejected. 

[0059] To extractthe background region 6, aseedfill algorithm is used that starts from the boundary of the collimation 
region as determined in the previous step. Appropriate seed points for B are found by considering a small band along 
each of the collimator edges and retaining all pixels whose intensity is smaller than the mean of the band. This approach 
avoids choosing pixels that belong to the diagnostic region as candidate seed pixels. The background region is then 

55 grown by considering all neighbouring pixels n,,/=1,...8 of each pixel pEBand adding g,to B if the intensity difference 
between p and q is smaller than some specified threshold. 
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Maximum likelihood estimation 

[0060] We simplify (13). by leaving out the multiplicative component s„of the image degradation process 



U* = argmin- X log e />(w,|/i,CT) 



[16) 



= arg min- £ „ log e p(U(x,y^,<?) 



B,U x,yzB 



= arg min- I Jog e p(N(x,y) - S A (x, y^n, <r) 



[0061] This equation is optimized by iteratively estimating the background B and finding parameters u., o and the 
components s A after differentation and substitution of p(u,\p.,a) by the Gaussian distribution (12). To find the solution 
for the multiplicative component, the same approach can be followed after logarithmic transforming the intensity values. 
[0062] The initial estimate for the background B is taken from the segmentation algorithm described higher. All other 
estimates for B are computed using a histogram based threshold algorithm. The threshold is defined as the smallest 
value of e satisfying 



where p^(n) is the probability that a point in image Op has value n and p., o are the mean and variance of the 
corrected pixels belonging to the previous background estimate. 

[0063] The maximum likelihood estimates for the parameters u. and a of 7, can be found by minimization of -Z,log e p 
(u ( l(j.,a). The expression for u. is given by the condition that 




(17) 




[0064] Differentiating and substituting p(u,lu.,a) by the Gaussian distribution (12) yields: 
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v = y,— = y Ar( *' ' - ) - Sa ^ ' yi } 



whereby,- is the spatial position of pixel / and n is the number of background pixels. The same approach can be 
followed to derive the expression for a : 



ff 2 = y (Pt^Jfr 



, (N(x h y;)-^ A (xi,yi)) 2 



[0065] Suppose that the spatially smoothly varying component s A can be modeled by a linear combination of K 
polynomial basis functions ty-^y,) 



u t = N(x t ,>>,)- £ Cj(P j (x t , y; ) 
j=l,...,K 



30 the partial derivative for c k of (1 6) set to zero yields 



Xj N(x h y { ) - fi -Z, c j<>j (wo = 0 v *- 



[0066] Solving this equation for {c^} does not seem very tractable, but combining all equations for all ft and introducing 
matrix notation simplifies the problem considerably 



where A represents the geometry of the image formation model, each of its rows evaluating one basis function 
§ k at all coordinates and R represents the residue image, i.e. the difference between the acquired image and the 
55 estimated background mean. In full matrix notation, the equation is 
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C = 



tf> 2 0c,)<M* 2 )k(* 3 )~ 




«, -11 



n 2 ~n 



where n,- is the intensity value of the acquired image at pixel (Xj.yj). Equation (18) is a least squares fit to the 
residue image. As least squares fit are sensitive to outliers, only entries in R which satisfy |n ( -^,|<2.5a are included to 
solve (18). 

Entropy Minimization 

[0067] Suppose that the image degradation components s A and s M can be modeled by a linear combination of K 
polynomial basis functions §j m - a (x,y) 



[0069] The optimal additive parameters a* and multiplicative parameters m are found by Powell's multidimensional 
directional set method and Brent's one-dimensional optimization algorithm (W.H. Press, S.A. Teukolsky, W.T. Vetterling, 
and B.P. Flannery. Numerical Recipes in C. Cambridge University Press, 1992.) 

[0070] The set of probabilities p(n) in (14) can be obtained by normalization of its histogram. In order to reduce the 
influence of random effects induced by discretizing the histogram, we use partial intensity interpolation at histogram 
formation. When transforming the image, an integer intensity value g is transformed to a real value g'. which in general 
lies between two integer values k and fc+1 . The histogram entries h(k) and h(/c+1) are updated by k+'\-g' and g'-k 
respectively. To obtain a smoother decline to the absolute minimum and to minimize the effects of local minima, the 
obtained histogram is blurred to: 



j=i k„ 



[0068] Equation (15) is reformulated as 
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where the parameter f was set to 3. 
Image Formation Models 

5 [0071] We have tested different image formation models which are summarized in Table 1 . The polynomial models 
are defined as 

$ v =c Q +c^x+c 2 y+c 3 x 2 +c 4 xy+c 5 y 2 +...+c y 



[0072] Models E,-,/=1 ,2 are included for the maximum likelihood estimation, model E 3 is the general image formation 
model while model Z 4 is derived from (2). Model S 5 is an approximation of model Z 4 where the different model param- 
eters are substituted with real values and higher orders are discarded where appropiate. Model E s is included for 
15 resemblance with model Z 2 ■ 



Table 1. Image formation models used for correction of 
inhomogeneities . 0' is a polynomial model of order v. 





Model 


Correction 


25 


E 


U = Nxj> v 




E 


U = N + <I) V 


30 


E 






e 


Jo? 

U = Af.exp(^V) 

<t>2~ 


35 








e 












E 





45 



Fourth embodiment 



so [0073] In a fourth embodiment according to the present invention, a statistical mixture model of the image is generated 
based on a plurality of K image regions. 

[0074] Each of these regions or classes may physically correspond to e.g. bone, soft tissue and direct exposure area. 
[0075] In the assumption of a normal mixture model, each class is represented by three unknown parameters: the 
proportion n k of image pixels, the mean value \i k and the variance a 2 . 
55 [0076] The set y collectively comprising all unknown parameters becomes: 
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V = 1t K , ,.. ., flR > °f >- . °K } 



[0077] The subset of parameters pertaining to class k is denoted as 



[0078] The image intensity histogram, denoting the probability distribution that a pixel / has intensity y,- is therefore 
a Gaussian mixture model 



/(y,lv)=S *Jtly,Wk ) 



The basic EM algorUhm 

[0079] The classical analytical method to estimate the parameter x|/ is to maximise the log-likelihood function for each 
of the parameters to estimate. 

[0080] The maximum likelihood estimates of each parameter can be solved from a system of equations which is 
non-linear in general and hence requires methods such as Newton-Raphson algorithm. 

[0081] The Expectation-Maximisation (EM) algorithm estimates the parameters y by adding segmentation labels z, 
(i represents pixel i and z, has a value k, k=1 ... K) , (so called non-observable data) to each of the grey values y, of 
the pixels (so called observable data). 

[0082] In each iteration of the EM algorithm the expectation step (E-step) estimates a segmentation label k\o each 
pixel / on the basis of parameter values y from the previous iteration and in the maximisation step (M-step) new pa- 
rameter values \|/ are computed on the basis of maximum likelihood, given the new segmentation labels associated 
with each of the newly assigned segmentation labels. 

The extended EM algorithm 

[0083] In the context of the present invention two modifications have been added to the EM algorithm to make it 
correcting for a bias field caused by global inhomogeneities in the imaging chain and to discard outliers due to local 
inhomogeneities. 

[0084] The global inhomogeneities in the image distort the assumed normal distribution of the pixel classes. 
[0085] Every pixel segmentation class is modelled as a normal distribution of which a sum of spatially correlated 
continuous basis functions is subtracted. 

[0086] Examples of such basis functions are orthogonal polynomials. Other orthogonal continuous functions may be 
used as well. 

[0087] The coefficients of the basis polynomials are added to the parameter set y which must be estimated 
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^{c v ,... > 7t K ,^,...,^ K ,(jl,...,al,c l ,...,c R } 

with the probability distribution for the pixels belonging to segmentation class k 



My\y*c)=-=L 

J27CG, 



with § r a Nx1 vector holding the polynomial function evaluation for the r-th basis polynomial at pixel location i 
(i=1...N). 

20 [0088] A further correction to the basic EM algorithm is to make it robust against outliers in the observed distribution 
of a segmentation class, caused by the presence of local defects (dust, scratch, pixel drop out...) in the recording 
member, said defects not being attributable to the global inhomogeneities. 

[0089] To this purpose each pixel class k is divided in a Gaussian class (which is distributed by the inhomogeneity 
and which is corrected by the bias function) and a rejection class. This rejection class is assumed to have a uniform 
25 distribution with probability density 5^ and contains a proportion ee[0,1] of the pixels. The probability distribution of 
pixel class k is therefore 

30 

Summary of [he extended EM algorithm 

[0090] The extended EM algorithm is summarised by the following formulas valid for iteration m: 
as E-slep: 

[0091] For each pixel class k, /c=/...Kand each pixel ;', /=/... N , compute 

40 /*CMv'l w) )*t' ) 



— ! — e*p(4« 2 ) 
2 



(m+1)_ 
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y, denoting the intensity values of pixel / 

\\i (m) Xhe set of statistical parameter describing class /rat iteration m 
n^the proportion of pixels in the image belonging to class k at iteration m 

f k the probability density function of intensity of pixels of class k denoting the conditional probability that pixel i 
has gray value y- t given parameters of class k 



a ' 'the variance of intensity of pixels belonging to class k at iteration m, 
k, a threshold on the Mahalanobis distance defined as 



tne p ro b a bjnty of pixels of class k being outliers, 
the probability of pixels inside class klo belong to the non-rejected group (i.e. not being an outlier). Because 
X k ± 0, this probability may be less than one, and hence 



[0092] At this stage, a segmentation of the image could be obtained by a hard classification, i.e. each pixel /' is 
assigned class frfor which p m+1> \s maximal, i.e. class pixel 



-tar'}- 



45 In the sequel of the EM algorithm, soft classification labels p^ +1> E[0...1] are used. 
M-step 

[0093] For each class k=1...K and for each coefficient c r r=1...R applied to the corresponding polynomial basis 
so function , compute 



55 
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N J m+l) 



.(«+!) 
"1 

.<«+!) 



= ( ^r r ( W)+l )^ ) -^r^(m + i) i?(OT+ i) 
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y\ - y\ 
y 2 -yt 



Z* 2 <m+l) Pk 



wherein 

|j/ m+ ^denotes the mean intensity value of pixels belonging to class kaX iteration (m+1), 

^2(m+i) denotes tne var j ance 0 f intensity value of pixels belonging to class k at iteration (m+1), after having 

corrected for the estimate of the bias field, 

C<™" 1 ) is a vector containing coefficients c p r=1 ...R applied to the corresponding polynomial basis function, 
A(i,i)=$ lr is the evaluation of the M-th polynomial basis function at pixel location ;' (matrix A thus represents the 

geometry of the bias field model), 

is a diagonal matrix of weights w/ m+1> , M...A/, 

with 

vJ m+1 ) the weight applied at pixel i in iteration (m+1). Said weight is the sum of the inverse of variance overall 
classes k, /c=1 ...K, each weighted with the probability of that class which is P^" +^ '' ^ |™ + '' , ■ 

is a residu image, the residu being the difference between the original image matrix yj 1=1. ..N and the 
corrected image matrix?'™*" at iteration (m+1). 

[0094] The equations of the extended EM algorithm reduce to the basic EM algorithm when no bias correction is 
performed (all Cf=0) or no outliers are taken into account (all X^and hence all t ik =1). 



Initialisation 



[0095] In order to start the iterations of the EM algorithm, an initial estimate \|/(0) for the parameter set y is needed. 
35 [0096] This is achieved by assigning each pixel /, i=1...N, to one of the classes k=1...K on the basis of intensity 
histogram slicing. 

[0097] This assignment involves the computation of , which is a hard assignment of probability 1 to one of the 
k possible class probabilities at pixel / and putting all other probabilities to zero. 
[0098] Furthermore, no outliers are assumed during initialisation, i.e. 1^=1 for all /. 
40 [0099] Therefor, the M-step in which the values \\r are computed can be executed immediately after initialisation. 
[0100] Therefore the initialisation step for which the iteration value m=0 does not present a true iteration step in the 
EM algorithm. 

[0101] To slice the histogram into K distinct initial pixel classes k=1...K, prior art techniques are used. In the context 
of the present invention, the histogram is smoothed and approximated with a higher order polynomial, after which the 
45 two or three most important maxima are determined. The intensity thresholds separating intensities of different classes 
are then determined as the intensities corresponding to the minima between these different maxima. 



Claims 



1. A method for correction of intensity inhomogeneities in a digital radiation image comprising (1) exposing an object 
to radiation emitted by a source of radiation, (2) recording a diagnostic radiation image of said object on a radiation- 
sensitive recording member, (3) reading the diagnostic radiation image that has been recorded on said recording 
member and converting the read image into a digital image representation, (4) generating a bias field comprising 

55 correction data, (5) correcting said digital image representation by means of said correction data, characterised 

in that said bias field is deduced from said digital image representation of said diagnostic radiation image. 

2. A method according to claim 1 wherein said digital image representation at least comprises a digital representation 
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of a direct exposure area and of a diagnostic area and (1) said digital image representation is subjected to image 
segmentation in order to extract data constituting an estimation of the direct exposure area and an estimation of 
the diagnostic areas and (2) said set of correction data is deduced from data constituting said estimation of the 
direct exposure area. 



3. A method according to claim 1 wherein said digital image representation at least comprises a digital representation 
of a direct exposure area and of diagnostic areas and (1) a mathematical model of a cause of said intensity inho- 
mogeneities is generated, (2) said digital image representation is subjected to image segmentation in order to 
extract data representing an estimation of the direct exposure area. (3) parameters of said mathematical model 
are deduced from image data corresponding to said direct exposure area, (4) a bias field is generated on the basis 
of the deduced parameters, (5) said image is corrected by means of said bias field so as to obtain a corrected 
image, (6) the corrected image is subjected to a stopping criterion, steps (2) to (6) are repeated on the corrected 
image until the stopping criterion is met. 



15 4. A method according to claim 3 wherein said cause of inhomogeneities is the Heel effect introduced by an X-ray 
source. 



5. A method according to claim 3 wherein said cause of inhomogeneities is the non-uniform sensitivity of a radiation 
detector. 



6. A method according to claim 3 wherein said stopping criterion is met when no significant changes occur in said 
estimation of the direct exposure area or in the parameters of said mathematical model. 

7. A method according to claim 1 wherein (1 ) a statistical model of said image is generated on the basis of intensity 
25 and spatial statistics of image regions in said image, (2) said digital image representation is subjected to image 

segmentation in order to extract data constituting an estimation of said image regions, (3) parameters of said model 
are estimated by means of data pertaining to said direct exposure area, (4) a bias field is generated, (5) said image 
is corrected by means of said bias field, (6) the result of step (5) is evaluated relative to a stopping criterion and 
steps (2) to (6) are repeated until said stopping criterion is met. 

30 

8. A method according to claim 7 wherein said image regions are confined to direct exposure areas. 

9. A method according to claim 7 wherein the method steps are applied separately to each of a number of image 
region classes jointly constituting the intensity histogram of the radiation image. 

35 

10. A method according to claim 7 wherein said statistical model is a Gaussian distribution of the direct exposure area 
and said parameters of the statistical model are the average value u. of the Gaussian distribution and the standard 
deviation a. 



40 11. A method according to claim 7 wherein said stopping criterion is reached when no significant changes occur in 
said estimation of image regions and/or no significant changes occur in the parameters defining said statistical 
model. 



12. A method according to claim 1 wherein said digital image representation at least comprises a digital representation 
of a direct exposure area and of diagnostic areas and (1) an information theoretic model of said image based on 
Shannon-Wiener entropy is generated, (2) said digital image representation is subjected to image segmentation 
in order to extract data representing an estimation of the direct exposure area, (3) the entropy in said model is 
extracted by means of data pertaining to said image regions, (4) a bias field is generated, (5) said image is corrected 
by means of said bias field. (6) the result of step (5) is evaluated relative to a stopping criterion and steps (2) to 
(6) are repeated until said stopping criterion is met. 

13. A method according to claim 12 wherein said stopping criterion is met when said entropy is minimal and no sig- 
nificant changes occur in the entropy value. 



55 
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